Method and apparatus for reducing multipath signal error using deconvolution

ABSTRACT

A deconvolution approach to adaptive signal processing has been applied to the elimination of signal multipath errors as embodied in one preferred embodiment in a global positioning system receiver. The method and receiver of the present invention estimates then compensates for multipath effects in a comprehensive manner. Application of deconvolution, along with other adaptive identification and estimation techniques, results in completely novel GPS receiver architecture.

ORIGIN OF THE INVENTION

The invention described herein was made in the performance of work under a NASA contract, and is subject to the provisions of Public Law 96-517 (35 U.S.C. § 202) in which the Contractor has elected not to retain title.

TECHNICAL FIELD

The present invention relates to signal error reduction, and more specifically to a method and apparatus for reducing Global Positioning Satellites (GPS) multipath signal error using deconvolution to determine the attitude and position of kinetic platforms.

BACKGROUND ART

The superiority of GPS receiver technology over other competing technologies, in terms of the cost, reliability, mass, size, and power considerations, as proven in its various applications such as navigation, spacecraft orbit determination, and surveying, has constantly prompted the extension of GPS to other important areas that require increasingly more exacting performance from the GPS system. One such area, for example, is the application of GPS to the attitude determination of aircraft and earth orbiting satellites.

For the GPS receivers to provide the requisite precision, the existing error sources in the GPS system must be eliminated or greatly reduced. At the current state of GPS technology, the most significant error source is the signal multipath propagation. For example, in spacecraft operational environments, the GPS signal is reflected from various structural components of the spacecraft and these reflected components are received by the GPS receiver along with the desired direct line-of-sight (LOS) path. The reflected signals differ from the desired LOS path signal in terms of their delays, amplitudes, and phases. The carrier phase tracking loop provides no inherent discrimination against the multipath signals and thus tracks the phase of the composite signal corrupted by multipath components. The resulting differential carrier phase estimation error can be orders of magnitude higher compared to the case of no multipath propagation in many GPS applications. For example,, the measurements obtained by the RADCAL satellite GPS-ADS (Attitude Determination System) experiments have shown that the differential range error in such environment is of the order of 1 cm corresponding to an attitude determination error of about 0.5 degree. Thus for the GPS receivers to provide precision pointing knowledge (order of 1 arcmin or better with 1 meter antenna baseline) or a differential range accuracy of about 0.3 mm or better, the multipath effects must be suppressed by orders of magnitudes. Similar accuracy may also be desirable in other GPS precision applications in the presence of multipath signals such as GPS based geophysical measurements and precision surveying.

Among the past approaches to deal with the multipath problem, one approach involves reducing the early-late delay spacing among the correlators in the GPS receiver code lock loop. However while this approach reduces the code range errors to some extent, it does not aid in the carrier phase measurements accuracy that is the basis of most GPS precision applications. Moreover even the reduction in the code range error is limited and if the early-late spacing is smaller than the initial delay error due to multipath signals (easily the case with many multipath situations), then the loop error can be very high and the loop may not even track.

In the work of VanNee et al., a set of nonlinear implicit equations are derived on the basis of maximum likelihood estimation theory. These equations contain the amplitude,, phase and delay of multipath signals. The paper proposes to solve these highly nonlinear equations in a recursive manner, but there is no explicit solution presented. The method has been implemented in a Nov Atel receiver using multiple correlators. Although this approach is useful for code phase measurement and a benign multipath environment, it may not be suitable for severe multipath environment and for precision GPS applications based on carrier phase measurements.

STATEMENT OF THE INVENTION

The present invention provides a method and apparatus for simultaneous estimation and compensation of the multipath estimation errors in both the carrier phase lock and the code delay lock loops in GPS receivers. The preceding is accomplished by disclosing new techniques for dealing with the multipath problem and includes these techniques in an improved GPS receiver. These techniques are based on the application of the optimal deconvolution approach in a somewhat unconventional manner, as compared to its application in other fields such as seismology and telecommunications. The proposed method consists of first estimating the impulse response of the effective multipath channel by a least squares algorithm. This step is followed by obtaining an inverse filter which equalizes the multipath channel response to the desired ideal multipath free response to the maximum extent possible within the specified constraints of the implementation complexity. From the equalized response, estimates of the true carrier phase and code delay can be made. The simulation results demonstrate that the proposed method is capable of reducing the multipath distortion by several orders of magnitude in an environment that is more severe than may possibly take place in any realistic precision GPS applications environment.

The invention also is embodied in a receiver comprising a series of correlators which generate demodulated signals from known reference signals, and the demodulated reference signals are then used to develop a discriminator function from which multipath error can be approximately determined. The present invention then compensates for the multipath error, and produces a final differential phase and code range estimate after compensating for the multipath error.

The present invention presents an optimum solution to the inherent tradeoffs among hardware/software implementation complexity, the extent of the multipath expected in the application, and the degree of multipath cancellation.

BRIEF DESCRIPTION OF THE DRAWINGS

The objects and features of the present invention, which are believed to be novel, are set forth with particularity in the appended claims. The present invention, both as to its organization and manner of operation, together with further objects and advantages, may best be understood by reference to the following description, taken in connection with the accompanying drawings.

FIG. 1 is a block diagram of the elements of a receiver in a preferred embodiment of the present invention;

FIG. 2 is a flow diagram of the channel impulse response estimator and equalizer of FIG. 1;

FIG. 3 is a graph illustrating the tracking of the Discriminator Function amplitude vs. delay in simulation example 1;

FIG. 4 is a graph illustrating the Discriminator response vs. delay in simulation example 1;

FIG. 5 is a graph of the normalized delay error magnitude vs. iterations for simulation example 1;

FIG. 6 is a graph of the pointing error vs. iterations for simulation example 1;

FIG. 7 is a graph of the pointing error vs. iterations for simulation example 1 with Delay>0.9;

FIG. 8 is a graph of the normalized delay error magnitude vs. iterations for simulation example 2;

FIG. 9 is a graph of the pointing error vs. iterations for simulation example 2;

FIG. 10 is a graph of the pointing error vs. iterations for simulation example 3; and

FIG. 11 is a graph of the normalized delay error vs. iteration number for simulation example 3.

DETAILED DESCRIPTION OF THE INVENTION

The following description is provided to enable any person skilled in the art to make and use the invention and sets forth the best modes contemplated by the inventors of carrying out their invention. Various modifications, however, will remain readily apparent to those skilled in the art, since the generic principles of the present invention have been defined herein specifically to provide a method and apparatus for reducing multipath error especially in GPS applications.

In the more conventional telecommunication applications of the deconvolution approach or the equalizer theory, the multipath propagation channel is modeled as, ##EQU1## where {u_(j) } represents the transmitted symbol sequence, {y_(j) } is the channel output sequence, and {h_(-q).sbsb.1, h_(-q).sbsb.1₊₁, . . . , h₀,h₁, . . . h_(q).sbsb.2 } represents the discrete channel impulse response. The additive noise sequence {v_(j) } is usually assumed to be a zero-mean white Gaussian. The basic process of deconvolution involves the estimation of the transmitted input sequence {u_(k) } on the basis of noisy observations {y_(k) } assuming that the discrete channel impulse response {h_(k) } is known to the receiver. In the case of unknown channel response, adaptive equalization techniques may be used wherein first an approximate estimate of {h_(k) } is obtained on the basis of a training sequence known in advance to the receiver and the channel output {y_(k) } and subsequently the estimate of {h_(k) } is refined adaptively with u_(k) replaced by its estimated/detected version in the adaptive algorithm. From the real time estimate of {h_(k) }, a time-varying inverse filter is derived which then filters {y_(k) } to obtain the estimate of {u_(k) }. In practice the two steps of estimating {h_(k) } and then finding the corresponding inverse filter are combined into a single step of finding directly the adaptive equalizer coefficients. The problem wherein the adaptive equalization is achieved without any training sequence is relatively more difficult and has also received considerable attention in the literature.

In the following, the deconvolution approach is extended to the problem of multipath elimination in the GPS receiver code tracking and carrier phase lock loops. Referring to the drawing of GPS receiver multipath compensator 100 shown in FIG. 1, the antenna 1 receives the GPS signal which is filtered and amplified by the RF front end 2. In the absence of the multipath signals, the output 3 of the RF front end 2 is given by: ##EQU2## where the receiver noise is not considered in the first instance, A_(c) is the received signal amplitude, ω_(c) is the GPS signal frequency, and a(t)=±1 is the pseudo-random code waveform that phase modulates the carrier. It is assumed that the data modulation is removed in a decision-directed manner. In the presence of N multipaths in addition to the direct line-of-sight path, the input signal may be characterized as

    s.sub.m (t)=s(t)+α.sub.1 s(t-τ.sub.1)+. . . +α.sub.N s(t-τ.sub.N)                                          (3)

where α_(i) and τ_(i) denote respectively the amplitude and delay of the i th multipath for i=1,2, . . . ,N, and s(t) is the direct path or line-of-sight component of s_(m) (t). After substituting (2) in (3), the composite received signal may be expressed as follows: ##EQU3##

This signal s_(m) (t) is translated from the carrier frequency to an intermediate frequency (IF) by converter 4 to yield IF signal s_(IF) (t) at an intermediate frequency ω₁. Thus, equation (4) can then be represented as ##EQU4## In the conventional delay-lock discriminator the signal in (5) is correlated with the reference signals ##EQU5## where t is the delay tracking error, τ_(d) is the offset delay referred to as early and late correlator delay and θ_(c) is the reference oscillator phase. The resulting correlation functions are given by

    s.sub.IF (t).s.sub.L (t)=a.sub.c α(t).α(t-τ-τ.sub.d)cos(θ.sub.m.sbsb.0 -θ.sub.C)+. . . +A.sub.c α.sub.N α(t-τ.sub.N).α(t-τ-τ.sub.d).cos(θ.sub.m.sbsb.N -θ.sub.C).                                        (7)

The resulting signal denoted by R_(1i) (t) may be written in the form below.

    R.sub.1i (τ)=h.sub.0i R.sub.c (τ+τ.sub.d)+h.sub.1i R.sub.c (τ+τ.sub.d -τ.sub.1)+. . . +h.sub.Ni R.sub.c (τ+τ.sub.d -τ.sub.N)h.sub.ki =α.sub.k A.sub.c cos(θ.sub.m.sbsb.k -θ.sub.c);k=0,1, . . . ,N  (8)

In equation (8) R_(c) (τ) represents the code autocorrelation function and α₀ =1. Similarly the corresponding signal obtained with s_(L) replaced by S_(E) in (7) and denoted by R_(2i) (t) may be written in the following form.

    R.sub.2i (τ)=h.sub.0i R.sub.c (τ-τ.sub.d)+h.sub.1i R.sub.c (τ-τ.sub.d -τ.sub.1)+. . . +h.sub.Ni R.sub.c (τ-τ.sub.d -τ.sub.N)                          (9)

The difference between R_(1i) and R_(2i), the so called discriminator function D_(i) (τ) is now given by

    D.sub.i (τ)=h.sub.0i g.sub.c (τ)+h.sub.1i g.sub.c (τ-τ.sub.1)+. . . +h.sub.Ni g.sub.c (τ-τ.sub.N) (10)

where

    g.sub.c (τ)=R.sub.c (τ+τ.sub.)-R.sub.c (τ-τ.sub.d

is the discriminator function in the ideal case. Clearly if α_(k) =0 for all k≠0 then D_(i) (τ)=g_(c) (τ) corresponding to the case of no multipath propagation. The standard delay lock loop converges to the solution τ_(m) of the equation D_(i) (τ)=0 instead of converging to 0 which is the solution of g_(c) (τ)=0. The multipath error τ_(m) can be excessive depending upon the actual multipath environment encountered as will be illustrated by several simulation examples.

In the approach of the present invention the discriminator function is measured first and then with equation (10) is used to estimate all the unknown variables including the multipath delays, amplitudes, and phases and the estimates of the multipath errors both in the phase lock and code lock loops via deconvolution. The multipath errors are then compensated for in arriving at the final differential phase and code range estimates.

To accomplish these objectives, multiple correlators corresponding to multiple values of τ_(d) are used in the present invention instead of just two correlators. Returning to FIG. 1, signal 5 is correlated with carrier signal 43 C_(i) (t) in correlator 6, and the output 7 is correlated in correlators 8,12,16 with signals 54, 53, and 55, respectively. The outputs 9,13, and 17 are received by integrate/dump units 10, 14, and 18 and each output is integrated over T seconds to generate the correlator functions represented in equation (8). In the same manner, all the correlator functions R_(i) (τ+iΔ) for all integers i in the range of -Q₁ to Q₂ are generated using additional correlators and integrate/dump units. Output 11, R_(i) (τ+Q₂ Δ) is the same as R₁₁ (τ) of equation (8) with τ_(d) in (8) replaced by Q₂ and the carrier signal 43 C_(i) (t) is given by

    C.sub.i (t)=2cos (ω.sub.1 t+θ.sub.c).          (11)

The multipath resolution parameter Δ is selected sufficiently small to provide the required multipath resolution and the number of correlators (Q₁ +Q₂ +1) is dependent on the extent of multipath delay spread and Δ. The linear combiner 20 receives output 11, 15, and 19 from the respective integrate/dump units 10,14, and 18, as well as any other integrate/dump units (not shown) and generates the sampled values of the function D_(i) (τ) given in (10) by differencing pairs of its input, that is,

    D.sub.i (τ+kΔ)=R.sub.i (τ+(k+I.sub.d)Δ)-R.sub.i (τ+(k-I.sub.d)Δ);k=-M.sub.i, . . . , 0, . . . , M.sub.2. (12)

where τ_(d) is selected equal to I_(d) Δ for some integer I_(d) and M₁ and M₂ values are dictated by the extent of multipath. Thus, in terms of M₁ and M₂, the Q₁ and Q₂ are given by Q₁ =M₁ +I_(d), Q₂ =M₂ +I_(d). The output 36 from the linear combiner 20 D_(i) is the (M₁ +M₂ +1) vector whose k th component is given by D_(i) (τ+kΔ) in (12).

The signal 5, S_(IF) (t) is correlated in 21 with the carrier signal C_(q) (t) 45 which is the quadrature phase version of C_(i) (t) and is given by:

    C.sub.q (t)=-1 sin (ω.sub.1 t+θ.sub.c)         (13)

Correlating output 22 of the correlator 21 with signal a(t-τ) in correlator 27 and integrating the output 28 over an interval of T seconds in the integrate/dump unit 29 generates R_(q) (τ) which is related to the parameters of interest by equation (14) as

    R.sub.q (τ)=h.sub.0q R.sub.c (τ)+h.sub.1q R.sub.c (τ-τ.sub.1)+. . . +h.sub.Nq R.sub.c (τ-τ.sub.N) h.sub.kq =α.sub.k A.sub.c sin(θ.sub.m.sbsb.k -θ.sub.c);k=0,1, . . . ,N                                                      (14)

In the same manner all of the correlator functions R_(q) (τ+iΔ) for all integers i in the range of -Q₁ to Q₂ are generated using additional correlators exemplified by 23,31 and whose outputs 24,32 are received by integrate/dump units 25,33 as shown in FIG. 1. Outputs 26,30,34 are transferred to the linear combiner 35 which generates (M₁ +M₂ +1) vector D_(q) whose k th component D_(q) (τ+kΔ) is given by:

    D.sub.q (τ+kΔ)=R.sub.q (τ+(k+I.sub.d)Δ)-R.sub.q (τ+(k-I.sub.d)Δ; k=-M.sub.1, . . . , 0, . . . , M.sub.2. (15)

From equation (15) it follows that D_(q) (τ) may be expressed as:

    D.sub.q (τ)=h.sub.0q g.sub.c (τ)+h.sub.1q g.sub.c (τ-τ.sub.i)+. . . +h.sub.Nq g.sub.c (τ-τ.sub.N). (16)

For the purpose of estimating and filtering, the outputs 36 and 37 are combined into the vector D given by

    D=D.sub.i +jD.sub.q                                        (17)

where j=√-1 and the k th component of D denoted by D^(k) is given by

    D.sup.k =D(τ+kΔ)=D.sub.i (τ+kΔ)+jD.sub.q (τ+kΔ)                                          (18)

From equations (10), (16), and (18) it follows that

    D(τ)=h.sub.0 g.sub.c (τ)+h.sub.1 g.sub.c (τ-τ.sub.1)+. . . +h.sub.N g.sub.c (τ-τ.sub.N)+n(τ) h.sub.k =h.sub.ki +jh.sub.kq ; k=1, . . . , N.                                         (19)

Here, D(τk) is the complex-valued discriminator function, and n(t) represents the noise at the correlator output corresponding to the noise at the receiver input.

FIG. 2 illustrates a flow chart of the steps executed by the channel impulse response estimator and equalizer 38. Step 200 involves obtaining the two discriminator functions D_(i) and D_(q) from the linear combiners 20,35 and solves for possible channel impulse responses h^(k) based on delay kΔ. The output from step 200 is a series of channel impulse responses h^(-N) through h^(N), which are then processed in step 210. Step 210 selects the impulse response which best fits the measurement for the inputs D_(i) and D_(q), denoted as h^(k).sbsp.0. In step 220, the corresponding deconvolution filter coefficients ƒ are computed from the results of step 210, and the coefficients are used in step 230 to remove the multipath contribution by convolving measurements {D^(j) } with {ƒ_(j) }. The resultant discriminator function D_(eq) is then solved in step 240 for Re{D_(eq) (τ)}=0 by linear interpolation of the sampled version of D_(eq). The delay obtained, τ is used in step 250 to solve for the channel impulse response h.sup.ƒ and to compute the argument (phase) of the central coefficient h₀.sup.ƒ in the final step 260.

Turning to the preferred solution, in FIG. 1 estimator 38 processes the input vectors D_(i) and D_(q) to estimate the multipath channel impulse response, obtain the coefficients of a deconvolution filter, and from these derive an estimate of the delay τ and the phase difference θ=(θ_(m0) -θ_(C)) such that these estimates are not corrupted by the multipath signals. For an understanding of the estimator 38, Equation (19) is rewritten in its discrete form (20), that is, substitute τ=jΔ where j is an integer into equation (19): ##EQU6## Here y_(j) =D(jΔ); g_(j) =g_(c) (jΔ); n_(j) =n(jΔ); and (-M₁,M₂) represents the interval over which the measured discriminator function is significant. Note that equation (20) is of somewhat more general form than (19) in that this form may also include negative values of the multipath delays. The delays τ₁, . . . , τ_(N) are approximated by integer multiples of Δ which can be achieved by selecting Δ sufficiently small. A least squares estimation technique is now applied to estimate the channel impulse response coefficients {h_(k) } and the corresponding deconvolution or equalizer filter coefficients from {y_(j) }. Letting

    h.sup.T = h.sub.-q.sbsb.1, . . . , h.sub.-1, h.sub.0, h.sub.1, . . . , h.sub.q.sbsb.2 !

    g.sub.j.sup.H = g.sub.j+q.sbsb.1, . . . , g.sub.j+1, g.sub.j, g.sub.j-1, . . . , g.sub.j-q.sbsb.2 !; j=-M.sub.1, . . . ,0, . . . , M.sub.2

the least squares estimate of the unknown channel impulse response vector is given by ##EQU7## where σ_(n) ² =E n_(j) ² ! and I is the identity matrix. Now denoting by χ₀ the truncated (or zero padded depending upon whether q_(i) >K_(i) or vice versa; i=1,2) version of h, i.e., with

χ₀ = h_(K).sbsb.1, . . . , h₀, . . . h_(-K).sbsb.2 !

where K₁ and K₂ are the selected numbers of the deconvolution filter coefficients, and by χ_(j) the j times shifted version of χ₀,

X_(j) = h_(K).sbsb.1_(+j), . . . , h_(j), . . . , h_(-K).sbsb.2_(+j) !

the optimum parameter vector f is given by ##EQU8## Note that in the definition of χ_(j) the entries for h_(j) are set equal to zero if j falls outside the interval -q₁,q₂ !. Also the constant κ is given by ##EQU9## The equalized channel response z is obtained by convolving the two sequences {h_(-q).sbsb.1, . . . ,h₀, . . . ,h_(q).sbsb.2 } and {ƒ_(-K).sbsb.1, . . . ,ƒ₀, . . . ,ƒ_(K).sbsb.2 }. The elements of z are denoted by {Z₋(q.sbsb.1^(+K).sbsb.1.sub.), . . . ,z₀, . . . ,z.sub.(q.sbsb.2_(+K).sbsb.2.sub.) }. Therefore the equalized discriminator response obtained by convolving {y_(j) } with {ƒ_(j) } sequence is given by ##EQU10## For the case of perfect equalization z₀ =1, z_(i) =0 for i≠0 and D_(eq) (τ)=g_(c) (τ). In practice Re{D_(eq) (τ)}≈g_(c) (τ) and the solution of the equation

    Re{D.sub.eq (τ)}=0                                     (24)

is the estimate of the delay τ. When the initial delay uncertainty τ is smaller compared to Δ/2 in magnitude, the j th components of the vector D_(i) and D_(q) at the input of 38 are approximately equal to the required real and imaginary parts of y_(j) in equation (20), and thus the deconvolution of the multipath channel is achieved by replacing y_(j) by the j th component of the vector D.tbd.D_(i) +√-1D_(q).

Delay Estimation

When the initial delay uncertainty τ is not small compared with Δ, it is explicitly taken into account by the estimator 38. Denoting the initial uncertainty by τ.sub.ε -τ.sub.ρ, the equation (19) is modified by replacing τ by (τ.sub.ε -τ.sub.ρ), yielding

    D(τ)=h.sub.0 g.sub.c (τ.sub.ε -τ.sub.ρ)+h.sub.1 g.sub.c (τ.sub.ε -τ.sub.1 -τ.sub.ρ)+. . . +h.sub.N g.sub.c (τ.sub.ε -τ.sub.N -τ.sub.ρ)+n(τ). (25)

In (25) τ.sub.ρ represents the negative of the signed integer part of (τ/Δ) and τ.sub.ε is the fractional part of (τ/Δ) such that

    τ=τ.sub.ε -τ.sub.ρ =τ.sub.ε -k.sub.0 Δ; |τ.sub.ε |<Δ/2

for some signed integer k₀ and with D selected above. The equation (25) may be rewritten as

    D(τ)=h.sub.0 g.sub.k.sbsb.0 (τ.sub.ε)+h.sub.1 g.sub.k.sbsb.0 (τ.sub.ε -τ.sub.1)+. . . +h.sub.N g.sub.k.sbsb.0 ·(τ.sub.ε -τ.sub.N)+n(τ) g.sub.k.sbsb.0 (τ.sub.ε)= g.sub.c (τ.sub.ε -k.sub.0 Δ)                                         (26)

Assuming that |k₀ |<N for some integer N, then k₀ is the integer k that minimizes the following index ##EQU11## where h^(k) is the channel response obtained on the basis of replacing k₀ by k in (26),

    g.sup.k (τ)= g.sub.k (τ.sub.ε),g.sub.k (τ.sub.ε -τ.sub.1), . . . ,g.sub.k (τ.sub.ε -τ.sub.N)!.sup.T (28)

and g_(k) (τ.sub.ε) is approximated by g_(c) (-kΔ). In the estimation of h^(k) for any integer k, the discretization and estimation procedure of (20)-(22) is followed. Thus the estimator 38 minimizes the following index ##EQU12## where g₀ ^(k) is the vector of dimension (q₁ +q₂ +1) whose lth component is equal to g_(c) (-(l+k)Δ), g_(j) ^(k) is equal to the j time shifted version of g^(k), and h^(k) is the solution of equation (21) with g_(j) replaced g_(j) ^(k).

If k₀ is the solution of the minimization in (29), then the overall channel estimate is given by h=h^(k).sbsp.0. Substitution of h in (22) provides the coefficients of the inverse filter and the solution of (24) yields τ, signal 46 in FIG. 1. Note that if it is known that the multipath error τ is smaller than Δ/2 in magnitude, then this additional procedure is not necessary. However for k₀ >0, considerable error can otherwise ensue in the estimate of h.

Multipath Phase Estimation

With τ denoting the estimate of the multipath delay error, one now solves equation (20) for h but with g_(j) equal to g_(c) (jΔ-τ) by least squares estimation procedure. Denoting the resulting least squares solution by h.sup.ƒ, then the multipath phase error estimate is given by the argument of the zeroth component of the estimated impulse response, i.e., θ=+arg (h₀.sup.ƒ) which corresponds to signal 39 in FIG. 1. This error is compensated from the carrier phase-lock loop phase estimate in the following estimation cycle.

Output 39 of estimator 38 is input to a loop filter 40, which may be a digital filter similar to the ones used in phase-locked loops. In its simple form it is a multiplier by a gain (1/G_(c) T) where is G_(c) the sensitivity of the numerically controlled oscillator (NCO) 42. Output 41 of the carrier loop filter 40 controls the phase of the carrier signal 43 (C_(i) (t)) at the output of NCO 42. Signal 43 is inputted into phase shifter 44 which generates a π/2 phase shifted signal 45 (C_(q) (t)). Output 46 of estimator 38 is input to a code loop filter 47, which may be a digital filter similar to the ones used in delay lock loops or phase lock loops. In its simple form it is a multiplier by a gain (1/G_(d) T) where G_(d) is the sensitivity of the NCO 49. Output 48 of the code loop filter 47 controls the phase (delay) of the output 50 of the code NCO 49. The output of code NCO is inputted to a code generator 51 whose output is the reference code a(t-τ+Q₁ Δ) which is the binary pseudo random code and is the delayed version of the code a(t) appearing in the received signal 3 (S_(m) (t).

The delay τ of this code is controlled by signal 48 at the input of code NCO 49. The code delay signal 52 is inputted to a tapped delay line 59 implemented by a shift register (not shown). The (Q₁ +Q₂ +1) outputs of this tapped delay exemplified by 52,53, and 54 are input to various correlators such as 8, 12, and 16, and to the correlators 23, 27, and 31.

In the close-loop estimator, τ is the estimate of the delay error τ between the code phase of the direct component of the incoming signal 5 and the phase of locally generated signal 53 at one sampling time instance earlier. This error is corrected via the signal 48 at the input to the code NCO 49. With the recursive procedure in the limit the error τ will approach zero and the code signal 53 will achieve lock with the direct path component of the incoming signal, thus eliminating the effect of multipath signal components present in signal 5.

Similarly, θ is the estimate of the phase error between the phase of the direct path component of the incoming signal 5 and the local reference signal 43 (C_(i) (t)), and by generating signal 41 via lock loop filter 40, the control of the phase of carrier NCO 42 is such as to drive this error to zero. With the recursive procedure the carrier NCO output signal 43 will achieve phase lock with the direct path component of the incoming signal, thus eliminating the effect of the multipath signal components on the carrier phase tracking. The accumulated value of signal 41 multiplied by G_(c) T or the accumulated value of signal 39 when the carrier loop filter 40 is a multiplier with gain (1/G_(c) T), is the estimate of the difference between the phase of the incoming signal direct path component and the initial carrier NCO signal phase. The value of signal 46 when the code loop filter is a multiplier with gain (1/G_(d) T), is the estimate of the delay between the code phase of the incoming signal direct path component and the initial code phase of the code NCO output 50.

For the case of GPS signals, this delay difference multiplied by the velocity of light is equal to the multipath error free estimate of the code phase range. Similarly the carrier phase difference multiplied by (λ/2π) where λ is the carrier signal wavelength, is equal to the multipath error free estimate of the carrier phase range, also known as the accumulated delta range.

SIMULATION RESULTS

In this section some simulation examples are presented depicting the performance of the proposed multipath cancellation algorithm. First some additional notations are introduced to present these results. Let T and T_(c) denote the sampling period for the signal processing system and the code chip period respectively. Also let τ_(d) denote the normalized delay τ_(d) /T_(c) and g be the sample signal-to-noise power ratio given by γ=(P_(c) T/N₀) where P_(c) =A_(c) ² /2 is the signal power received by the direct or line-of-sight path and N₀ is the one-sided power spectral density of the receiver thermal noise. Further denote by a and θ_(m) the vectors consisting of the amplitudes and phases of the discrete channel response h in equation (26). One motivation behind this invention has been the GPS application to precision attitude determination of LEO (low earth orbit) spacecrafts. Therefore the achievable carrier phase estimation errors are interpreted in terms of the equivalent attitude pointing errors. It may be easily seen that for 1 meter antenna baseline the conversion factor between the carrier phase error in radians and the corresponding pointing error in arcmin is approximately equal to 109.4. Several simulation examples are presented below in terms of the notations introduced in this section. Recall that (q₁ +q₂) is equal to the number of discrete multipaths , (K₁ +K₂ +1) is the number of filter taps and g is the sample SNR.

Example 1

q₁ =0; q₂ =8; K₁ =12; K₂ =4; γ=10³ ;τ_(d) =0.5

a= 1 0.2 0.5 0.8 0.9 0.7 1 0.2 0.9!

θ_(m) = 8 0 -0.5 0.7 1.2 0.8 -0.6 1.3 0.2 -1.1!

FIG. 3 plots the ideal discriminator response g_(c) (τ), distorted response D(t) of (19), and the discrete version of the equalized response D_(eq) (τ). Clearly while the zero-crossing of D(t) is about 0.4 T_(c), the zero-crossings of both g_(c) (τ) and D_(eq) (τ) are equal to zero. Reducing the tap spacing to τ_(d) =0.1 does not make any significant difference in the zero-crossing of the function D(t) as shown in FIG. 4, which plots the various discriminator functions for the case of τ_(d) =0.1. FIG. 5 shows the convergence of the delay error as a function of the number of samples processed with an initial normalized delay of -0.498 chips.

It can be seen that the steady state error without multipath correction will remain approximately equal to 0.4 T_(c) while with the correction algorithm the error is approximately equal to 0.001 T_(c). Note also that with a 50 dB-Hz CNR at the receiver input this period of 100 iterations corresponds to only 1 sec of real time as per the definition of g. FIG. 6 shows the carrier phase error expressed in arcmin for the GPS attitude determination application. Note that the initial error of more than 30 arcmin is reduced to about 0.1 arcmin.

Example 2

q₁ =0; q₂ =25; K₁ =35; K₂ =4; γ=10⁴ ; τ_(d) =0.1

a= 1 0.8 0.5 0.3 0.9 0.2 0.5 0.3 0.7 0.1 0.05 0 0.08 0.1 0.03 0.02 0.1 0.07 0.05 0.04 0.1 0 0.05 0.08 0.01!

θ_(m) = 0 0.5 0.7 -0.7 1 0.3 0.5 -0.6 1.5 1 0.8 0.2 0.1 0 1 -0.5 -0.3 0.6 0.1 0.5 -0.8 -0.5 0.9 1.2 -0.6 1!

FIGS. 7 and 8 depict respectively the convergence of carrier phase and the code delay errors versus the number of samples. In this example the sample SNR is equal to 40 dB and thus corresponds to an update period of 0.1 sec for a 50 dB-Hz CNR at the GPS receiver input. Note that the multipaths with delay greater than 9 D are weak in that their relative amplitudes are less than 0.1 of the direct path which makes the identification of these multipaths more difficult. As may be inferred from these figures a residual pointing error of less than 1 arcmin and a normalized delay error of 0.005 chips is obtained in less than 50 iterations. Since the amplitudes of multipaths for delays exceeding 9 D are relatively small, one may consider ignoring these and thus select a filter order that is much smaller than q. FIG. 9 shows the resulting error performance when the filter order K₁ is selected to be 10. As is apparent from the figure the residual error, while much smaller compared to the initial error of about 30 arcmin, is much higher than that achieved with K₁ =35 filter case.

Example 3

q₁ 0; q₂ =35; K₁ =35; K₂ =4; γ=10³ ; τ_(d) =0.1

a= 1 0.8 0.5 0.3 0.9 0 0.2 0.5 0.3 0.7 0.1 0.5 0.6 0.2 0 1 0.7 0.5 0.2 0.3 0.5 0.7 1 0.4 0.8 0.3 0.7 1 0.7 0.8 0.6 1 0.5 0.3 0 0.7!

θ_(m) = 0 0.5 0.7 -0.7 1 0.3 0.5 -0.6 1.5 1 0.8 0.2 0.1 0 1 -0.5 -0.3 0.6 0.1 0.5 -0.8 -0.5 0.9 1.2 -0.6 1 0.7 -1 -0.7 0.8 0.6 0.5 -0.8 1 0.4 -0.9!

FIGS. 10 and 11 plot the residual carrier phase and code delay errors respectively, for g!=30 dB. As may be inferred from these figures, for this case of very severe multipath propagation, the pointing error is reduced from about 32 arcmin to only 4.5 arcmin while the residual delay error is about 0.01 chips.

It will be understood that the embodiment described herein are merely exemplary and that a person skilled in the art may make many variations and modifications without departing from the spirit and scope of the invention. For example, one such modification comprises sequential time sharing among a relatively small number of correlators to obtain correlation signals over a much larger range of delay offsets. This is achieved by sequentially offsetting the code NCO phase over a range of delay offsets and computing correlation signals for each such offsets instead of computing all the correlation signals in parallel as described above in the preferred embodiment. Such sequential method may be attractive if the phase and delay of the input signal vary relatively slowly with time. All such variations and modifications are intended to be included within the scope of the invention as defined in the appended claims. 

We claim:
 1. A method for evaluating and compensating for multipath signal error in code delay and carrier phase comprising the steps of:(a) receiving an input signal comprising a direct path signal and a plurality of multipath signals; (b) generating a first set of correlation signals by correlating the input signal with a known reference carrier signal and a plurality of reference code signals; (c) generating a real component of a discriminator function by subtracting pairs of said first set of correlation signals; (d) generating a quadrature phase reference carrier signal obtained by shifting the known reference carrier signal by ninety degrees (e) generating a second set of correlation signals by correlating the input signal with the quadrature phase reference carrier signal and the plurality of reference code signals; (f) generating an imaginary component of the discriminator function by subtracting pairs of said second set of correlation signals; (g) computing a carrier phase difference estimate between a phase of the known reference carrier signal and a carrier phase of the direct path component of the input signal, and also computing a code delay offset estimate between a selected reference code signal phase and a code phase of the direct path component of the input signal; (h) processing the carrier phase difference estimate using a carrier loop filter, and updating the phase of the known reference carrier signal to minimize the carrier phase difference; (i) processing the code delay offset using a code loop filter, and updating the selected reference code signal to minimize the code delay offset; (j) repeating steps (b) through (i) in real time such that said carrier phase difference and said code delay offset are maintained at values less than predetermined limits; and (k) outputting a signal corresponding to an estimated direct path component of said input signal.
 2. The method as recited in claim 1 further including the step of converting the received input signal from a carrier frequency to an intermediate frequency prior to generating said first set of correlation signals.
 3. The method as recited in claim 1 wherein the known reference signal is generated by a numerically controlled oscillator.
 4. The method as recited in claim 1 wherein the plurality of reference code signals are generated by a numerically controlled oscillator, a code generator, and a tapped delay line.
 5. The method as recited in claim 1 in which step (g) further comprises the steps of:(g1) generating a number of sampled multipath channel impulse response estimates corresponding to a number of uniformly spaced values of delay offsets using a least squares estimation algorithm; (g2) selecting one of the channel impulse response estimates based on minimizing the magnitude of a square of the difference between a complex valued sample of a discriminator function, whose real and imaginary components are respectively equal to said values of said first and second set of discriminator functions and sample values generated with the delay offset among said uniformly spaced values of delay offset and a corresponding channel impulse response estimate; (g3) generating a vector comprising the coefficients of an optimum deconvolution filter with said channel impulse response estimate and its shifted versions; (g4) convolving said complex valued sampled of said discriminator function with said coefficients of the deconvolution filter to generate a sampled version of an equalized discriminator function; (g5) estimating the code delay offset as the value for which the real part of the equalized discriminator function is equal to zero; (g6) generating a fixed multipath channel impulse response estimate based on said code delay offset estimate; and (g7) estimating the carrier phase difference as an argument of the complex valued sample with an index equal to zero of said final multipath channel impulse response estimate.
 6. A signal receiving apparatus which compensates for multipath signal error comprising:means for receiving an input signal comprising a direct path signal and a plurality of multipath signals, said direct path signal and said plurality of multipath signals each comprising a carrier component and modulated by a code component; means for processing said input signal in first and second sets of parallel correlators to generate first and second sets or correlation signals, said correlators comprising means for correlating said input signal with a reference carrier signal in the first set of correlators and with a ninety degree phase shifted signal with respect to said reference carrier signal in the second set of correlators, wherein each set of correlators correlate using delayed signals derived from a common reference code signal, said first and second set of correlators comprising at least three correlators; means for combining the correlation signals generated within each of said first and second set of correlators to generate values of a first and second discriminator functions, said values of the first and second discriminator function comprising real and imaginary components of a complex-valued discriminator function; means for generating estimates of a phase error between a phase of the reference carrier signal and a carrier phase of the direct path component of the input signal; means for generating estimates of a delay offset between a selected reference code signal phase and a code phase of the direct path component of said input signal; carrier loop filter means for compensating the phase of the reference carrier signal based on the estimate of the phase error to minimize said phase error; code loop filter means for compensating the phase of the reference code signal based on the estimate of the delay offset to minimize said delay offset; and means for generating the direct path signal component of said input signal from said phase error and said delay offset.
 7. The signal receiving apparatus as recited in claim 6 further comprising:a carrier numerically controlled oscillator, and a code numerically controlled oscillator,wherein said carrier loop filter means compensates said carrier numerically controlled oscillator to recursively converge the reference carrier phase with the direct path carrier phase, and said code loop filter means compensates said code numerically controlled oscillator to recursively converge the selected reference code signal phase with the direct path code phase. 